function [S_stat, K_vec]= Robust_stat_calc(data,theta,Farray,giveMoments,giveWeight,giveVarianceEstim)
%Calculates S statistic, as well as vector of composite K statistics corresponding to the
%linear combinations in Farray. giveMoments is used to pass the moment and
%derivate functions, while giveWeight is used to pass the desired weighting
%matrix and giveVarianceEstim passes a function used to calculate
%covariance matrices

moments=str2func(giveMoments);
VarianceEstimator=str2func(giveVarianceEstim);
W_func=str2func(giveWeight);

obs=size(data,1);

%Calculate moments for each observation
gmat=moments(data,theta,0);
g=mean(gmat)';

%Calculate derivatives for each observation
q_array=moments(data,theta,1);

Sigma_base=gmat;
for n=1:size(q_array,3)
    Sigma_base(:,end+1:end+size(q_array,2))=q_array(:,:,n);
end

%Calculate joint asymptotic variance matrix of moments and Jacobian
Sigma=VarianceEstimator(Sigma_base);
Sigma_gg=Sigma(1:size(gmat,2),1:size(gmat,2)); %Variance of moments
Sigma_tg=Sigma(size(gmat,2)+1:end,1:size(gmat,2)); %Covariance of derivatives with moments
Sigma_tt=Sigma(size(gmat,2)+1:end,size(gmat,2)+1:end); %Variance of derivaties    
Sigma_gg_inv=Sigma_gg^-1; %Calculate and store inverse of variance of moments

%Calculate orthogonalized Jacobian Matrix
q_stack=Sigma_base(:,size(gmat,2)+1:end);
q=mean(q_stack)';
D_vec=q-Sigma_tg*Sigma_gg^-1*g; %Calculate vectorized version of orthogonalized Jacobian
D=reshape(D_vec,size(gmat,2),[]); %reshapes the vectorized orthogonalized Jacobian to the orthogonalized Jacobian

%Calculate weighting matrix
W=W_func(data,theta);

%Calculate S statistic of Stock and Wright
S_stat=obs*g'*Sigma_gg_inv*g;

%Calculates K statistic for linear hypotheses specified by Farray
K_vec=zeros(size(Farray,3),1);  
for n=1:size(Farray,3)
    F=Farray(:,:,n);
    F(:,sum(abs(F),1)==0)=[];
    K=obs*g'*W*D*(D'*W*D)^-1*F*...
            (F'*(D'*W*D)^-1*D'*W*Sigma_gg*W*D*(D'*W*D)^-1*F)^-1*...
            F'*(D'*W*D)^-1*D'*W*g;
    K_vec(n,1)=K;
end
end